library(estimatr); library(broom); library(dplyr)
library(rio); library(data.table); library(ggplot2)
library(ivDiag)


gsz_north <- import("metadata/ltp1F.dta") %>% setDT
gsz_south <- import("metadata/ltp2F.dta") %>% setDT
north_samp <- gsz_north[is.na(regione) | regione != 20]
south_samp <- gsz_south[sud == 1 & regist != 20 & dummyroma != 1]


Y1 <- "totassoc_p"  # non-profit organizations
Y2 <- "sede_aido" # organ donation
D <- "libero_comune_allnord"
Z <- "bishopcity"
weights <- "population"
controls <- c('altitudine', 'escursione', 'costal', 'nearsea',
              'population', 'pop2', 'gini_land', 'gini_income')
controls2 <- c('altitudine', 'escursione', 'costal', 'nearsea',
              'population', 'pop2', 'gini_land', 'gini_income',
              'capoluogo')              

# reduced form formula
rf1_fml <- as.formula(paste(Y1, "~", Z, "+", paste(controls, collapse = "+")))
rf2_fml <- as.formula(paste(Y2, "~", Z, "+", paste(controls, collapse = "+")))


north_rf1 <- lm_robust(rf1_fml,
    data = north_samp[dummyroma == 0 & totassoc_p > 0 & !is.na(citta_etrusca_allnord)], weights = population,
    se_type = "HC1")

north_rf2 <- lm_robust(rf2_fml,
    weights = population,
    data = north_samp[dummyroma == 0 & !is.na(citta_etrusca_allnord)],
    se_type = "HC1")


rf1_fml2 <- as.formula(paste(Y1, "~", Z, "+", paste(controls2, collapse = "+")))
rf2_fml2 <- as.formula(paste(Y2, "~", Z, "+", paste(controls2, collapse = "+")))

south_rf1 <- lm_robust(rf1_fml2, data = south_samp, weights = south_samp$population, se_type = "HC1")
south_rf2 <- lm_robust(rf2_fml2, data = south_samp, weights = south_samp$population, se_type = "HC1")

###########################
## Table A1 - reduced form
###########################

summary(north_rf1)$coefficients["bishopcity", 1:2]
north_rf1$nobs
summary(north_rf2)$coefficients["bishopcity", 1:2]
north_rf2$nobs
summary(south_rf1)$coefficients["bishopcity", 1:2]
south_rf1$nobs
summary(south_rf2)$coefficients["bishopcity", 1:2]
south_rf2$nobs

###########################
## Figure A11
###########################

## ZFS test
df <- north_samp[dummyroma == 0 & totassoc_p > 0 & !is.na(citta_etrusca_allnord)]
prior1 <- south_rf1 %>% tidy() %>% filter(term == "bishopcity") %>% select(estimate) %>% pull()
prior2 <- south_rf2 %>% tidy() %>% filter(term == "bishopcity") %>% select(estimate) %>% pull()

# Nonprofit
(ltz_out1 <- ltz(data = df, Y = Y1, D = D, Z = Z,
    controls = controls, weights = weights,
    prior = c(prior1, sqrt(0.2))))

p1 <- plot_ltz(ltz_out1, xlim = c(-1, 7))
ggsave("graphs/FigureA11a.pdf", p1, width = 4, height = 6)


# Organ donation
ltz_out2 <- ltz(data = df, Y = Y2, D = D, Z = Z,
    controls = controls, weights = weights,
    prior = c(prior2, sqrt(0.2)))

p2 <- plot_ltz(ltz_out2, xlim = c(-3, 4))
ggsave("graphs/FigureA11b.pdf", p2, width = 4, height = 6)
